module statistics
    implicit none
    real(8) :: sumk_mc, sumk_cmfd, sumk2_mc, sumk2_cmfd
    real(8) :: avgk_mc, avgk_cmfd, stdk_mc, stdk_cmfd, stdk_mcc
    integer, private :: nacc(2)
contains

subroutine init_statistics
    sumk_mc = 0.
    sumk2_mc = 0.
    sumk_cmfd = 0.
    sumk2_cmfd = 0.
    avgk_mc = 0.
    avgk_cmfd = 0.
    stdk_mc = 0. 
    stdk_cmfd = 0. 
    nacc = 0;

end subroutine

subroutine st_avgk(k_mc)
    real(8) :: k_mc
    integer :: nc
    nacc(1) = nacc(1) + 1
    nc = nacc(1)
    sumk_mc = sumk_mc + k_mc
    sumk2_mc = sumk2_mc + k_mc**2
    avgk_mc = sumk_mc/nc
    
    if (nc.ne.1) then
        stdk_mc = sqrt((sumk2_mc-sumk_mc**2/nc)/nc/(nc-1))
        stdk_mcc = sqrt((sumk2_mc-sumk_mc**2/nc)/(nc-1))
    end if 
end subroutine


subroutine st_avgkcmfd(k_cmfd)
    real(8) :: k_cmfd
    integer :: nc
    nacc(2) = nacc(2) + 1
    nc = nacc(2)
    sumk_cmfd = sumk_cmfd + k_cmfd
    sumk2_cmfd = sumk2_cmfd + k_cmfd**2
    avgk_cmfd = sumk_cmfd/nc
    
    if (nc.ne.1) then
        stdk_cmfd = sqrt((sumk2_cmfd-sumk_cmfd**2/nc)/nc/(nc-1))
    end if 
end subroutine

!subroutine add_statistics(k_mc, k_cmfd)
!    real(8) :: k_mc, k_cmfd
!    
!    nacc = nacc + 1
!    sumk_mc = sumk_mc + k_mc
!    sumk2_mc = sumk2_mc + k_mc**2
!    avgk_mc = sumk_mc/nacc
!    sumk_cmfd = sumk_cmfd + k_cmfd
!    sumk2_cmfd = sumk2_cmfd + k_cmfd**2
!    avgk_cmfd = sumk_cmfd/nacc
!    
!    if (nacc.ne.1) then
!        stdk_mc = sqrt((sumk2_mc-sumk_mc**2/nacc)/nacc/(nacc-1))
!        stdk_mcc = sqrt((sumk2_mc-sumk_mc**2/nacc)/(nacc-1))
!        stdk_cmfd = sqrt((sumk2_cmfd-sumk_cmfd**2/nacc)/nacc/(nacc-1))
!    end if 
!end subroutine    

end module